% Implements a repeated measures ANOVA with the factors height (3 levels)
% and weight (2 levels)

% Needs Eric's data files "pureTotal27-1.mat" and "ehiSCORE.mat"

% read hand switch data from Eric's hand selection data files
dataLowNoWeight = [totalBiases{:,1}]';
dataMidNoWeight = [totalBiases{:,2}]';
dataHighNoWeight = [totalBiases{:,3}]';
dataLowWeight = [totalBiases_weighted{:,1}]';
dataMidWeight = [totalBiases_weighted{:,2}]';
dataHighWeight = [totalBiases_weighted{:,3}]';

gripStrengthLeft = [freqData{:,8}]';
gripStrengthRight = [freqData{:,9}]';
gripStrengthDiff = gripStrengthLeft-gripStrengthRight;

% create data table
dataTable = table(dataLowNoWeight,dataMidNoWeight,dataHighNoWeight,...
    dataLowWeight,dataMidWeight,dataHighWeight,...
    'VariableNames',{'switchLowNoWeight' 'switchMidNoWeight' 'switchHighNoWeight' 'switchLowWeight' 'switchMidWeight' 'switchHighWeight'})

% create factorial condition table
FactorLevels = dataset({'L';'M';'H';'L';'M';'H'},{'N';'N';'N';'Y';'Y';'Y'},'VarNames',{'Height','Weight'})

% specify rm model
rm = fitrm(dataTable,...
    'switchLowNoWeight,switchMidNoWeight,switchHighNoWeight,switchLowWeight,switchMidWeight,switchHighWeight~1','WithinDesign',FactorLevels)

disp('Estimated Marginal Means');
margmeans{1} = margmean(rm,{'Height'});
margmeans{1}
margmeans{2} = margmean(rm,{'Weight'});
margmeans{2}
margmeans{3} = margmean(rm,{'Height' 'Weight'});
margmeans{3}

% Calculate rANOVA results for main effects and interaction
[ranovaResults,a,c] =  ranova(rm,'WithinModel','Height*Weight');
ranovaResults

% do sphericity test for 3 level factor Height and Interaction
mauchlySphericityTests = mauchly(rm,c)

% post-hoc for weight, simply to get mean of the differences
if ranovaResults(5,5).pValue < .05
    WeightposthocLSD = multcompare(rm,'Weight','ComparisonType','lsd')
else
    WeightposthocLSD = [];
end;

% post-hocs for main effect of height, Greenhouse-Geisser corr if Mauchly < .05
if mauchlySphericityTests(2,4).pValue > .05
    if ranovaResults(3,5).pValue < .05
        HeightposthocLSD = multcompare(rm,'Height','ComparisonType','lsd')
    else
        HeightposthocLSD = [];
    end;
else
    if ranovaResults(3,6).pValueGG < .05
        HeightposthocLSD = multcompare(rm,'Height','ComparisonType','lsd')
    else
        HeightposthocLSD = [];
    end;
end;

% post hocs to specify interaction
% compare weight levels at each height separately
tLabels = {'Low','Mid','High'};

if mauchlySphericityTests(4,4).pValue > .05
    if ranovaResults(7,5).pValue < .05
        [~,p(1,1),~,stats] = ttest(dataLowNoWeight,dataLowWeight);
        t(1,1) = stats.tstat;
        df(1,1) = stats.df;
        [~,p(2,1),~,stats] = ttest(dataMidNoWeight,dataMidWeight);
        t(2,1) = stats.tstat;
        df(2,1) = stats.df;
        [~,p(3,1),~,stats] = ttest(dataHighNoWeight,dataHighWeight);
        t(3,1) = stats.tstat;
        df(3,1) = stats.df;
        posthocInteractWeightDiff = table(p,t,df,'rownames',tLabels)
    else
        posthocInteractWeightDiff = [];
    end;
else
    if ranovaResults(7,6).pValueGG < .05
        [~,p(1,1),~,stats] = ttest(dataLowNoWeight,dataLowWeight);
        t(1,1) = stats.tstat;
        df(1,1) = stats.df;
        [~,p(2,1),~,stats] = ttest(dataMidNoWeight,dataMidWeight);
        t(2,1) = stats.tstat;
        df(2,1) = stats.df;
        [~,p(3,1),~,stats] = ttest(dataHighNoWeight,dataHighWeight);
        t(3,1) = stats.tstat;
        df(3,1) = stats.df;
        posthocInteractWeightDiff = table(p,t,df,'rownames',tLabels)
    else
        posthocInteractWeightDiff = [];
    end;
end;

save 'handSelection_Biases_rANOVA' dataTable margmeans rm ranovaResults mauchlySphericityTests WeightposthocLSD HeightposthocLSD posthocInteractWeightDiff;